# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Collective victimhood beliefs and conflict-related attitudes: A meta-analysis
# Marko Kljajic, Nadav Shelef, and Ethan vanderWilden

# Last Replicated: August 6, 2025

# Helper functions
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


# A. Function that plots the raw data (+ 95% CIs)
meta_plot <- function(data, range = 3.2, jump = 1, sepColor = F, includeEst = F){
  
  #Separate into inclusive and non-inclusive CVB estimates
  data$inc <- ifelse(data$treatment == "Inclusive victimhood", "Inclusive", "Non-Inclusive")
  
  #Extract overall estimate
  numResults <- meta3L(y = data$d.directional, v = data$var.d, cluster = data$article.id)
  myEst <- round(summary(numResults)$coefficients[1,1], 2)
  
  #Plot results
  myplot <- ggplot(data = data) +
    geom_pointrange(aes(x = d.directional, y = reorder(estimate.id, d.directional), 
                        xmin = d.directional - 1.645*sd.d, xmax = d.directional + 1.645*sd.d, color = inc),
                    alpha = 0.4, shape = 15) +
    geom_vline(xintercept = 0, linetype = "dashed")+
    labs(x = "Effect size (90% CI)", y = "Estimate") +
    scale_x_continuous(limits = c(-range, range), breaks = seq(-round(range), round(range), jump)) +
    theme_bw() +
    theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
          panel.grid.major.y = element_blank(), axis.title = element_text(size = 14)) +
    scale_color_manual(values = c("black", "black"))
  
  #separate by inclusive and non-inclusive victimhood if we want color coding
  if(sepColor == F){
    myplot <- myplot + theme(legend.position = "none")
  } else{
    myplot <- myplot + scale_color_manual(values = c("darkblue", "darkred"))+
      theme(legend.position = "top", legend.title = element_blank(), legend.text = element_text(size = 14))
  }
  
  #if we want to add the overall meta-estimate, add corresponding vertical line
  if(includeEst == T){ myplot <- myplot + geom_vline(xintercept = myEst, color = "black")}
  
  
  return(myplot)
}

# B. Function that returns results table (6 x 5 or 6 x 5) from 3-level meta analysis (with combinations)
fill_metatable3 <- function(data, exclude_mechs = F){
  
  #### A. create dataframe to store results
  if (exclude_mechs == F){ # if we are including emotions and cognitive frames...
    estimates <- as.data.frame(matrix(nrow = 5, ncol = 7))
    row.names(estimates) <- c("All victimhood", "General", "Competitive", "Exclusive", "Inclusive")
    names(estimates) <- c("All outcomes","Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment", "Perspectives", "Emotions")
    
    #specify the victimhood types and outcomes
    victimhood <- c("Collective victimhood", "Competitive victimhood", "Exclusive victimhood", "Inclusive victimhood")
    outcomes <- c("Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment", "Cognitive Frame", "Emotions")
  } else {
    estimates <- as.data.frame(matrix(nrow = 5, ncol = 5))
    row.names(estimates) <- c("All victimhood", "General", "Competitive", "Exclusive", "Inclusive")
    names(estimates) <- c("All outcomes","Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment")
    
    #specify the victimhood types and outcomes
    victimhood <- c("Collective victimhood", "Competitive victimhood", "Exclusive victimhood", "Inclusive victimhood")
    outcomes <- c("Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment")
  }
  
  
  #### B. Make table 4 x 6 (5x7 with totals) table [or 4 x 4 (5x5 with totals)]
  #         - logic: run through each victimhood and outcome combination. 
  #                  Subset data to just those estimates. Run 3-level meta analysis and store result
  for (i in 1:length(victimhood)){
    for (j in 1:length(outcomes)){
      #subset by victimhood category
      subdata <- subset(data, data$treatment == victimhood[i])
      
      #subset by outcome category
      subdata <- subset(subdata, subdata$outcome == outcomes[j])
      
      #if there are no estimates, then fill the table in as an NA
      if (nrow(subdata) == 0){
        estimates[i+1,j+1] <- NA
      } else if (nrow(subdata) == 1){
        estimates[i+1,j+1] <- paste(round(subdata$d.directional, 2), " (", round(subdata$sd.d, 2), ")", sep = "")
      } else {
        metaresult <- meta3L(y = subdata$d.directional, v = subdata$var.d, cluster = subdata$article.id)
        metaresult <- summary(metaresult)$coefficients
        
        #if metaresult cannot produce SEs because of failed openMX output, take SE from standard REML
        if(is.na(metaresult[1, 2]) == T){
          retry <- rma(yi =subdata$d.directional, vi = subdata$var.d)
          metaresult[1,2] <- retry$se
          metaresult[1,6] <- retry$pval
        }
        
        #store estimate, se, and significance in the estimates table
        estimates[i+1,j+1] <- paste(format(round(metaresult[1,1], digits = 3), nsmall = 3), 
                                    " (", format(round(metaresult[1,2], digits = 3), nsmall = 3), ")", 
                                    ifelse(metaresult[1,6] < 0.05, "*", ""), sep = "") 
      }
    }
  }
  
  
  #Fill in the 'TOTAL' column for each collective victimhood belief
  for (i in 1:length(victimhood)){
    subdata <- subset(data, data$treatment == victimhood[i])
    
    #if there are no estimates, then fill the table in as an NA
    if (nrow(subdata) == 0){
      estimates[i+1,1] <- NA
    } else if (nrow(subdata) == 1){
      estimates[i+1,1] <- paste(round(subdata$d.directional, 2), " (", round(subdata$sd.d, 2), ")", sep = "")
    } else{
      metaresult <- meta3L(y=subdata$d.directional, v=subdata$var.d, cluster=subdata$article.id)
      metaresult <- summary(metaresult)$coefficients
      
      #store estimate, se, and significance in the estimates table
      estimates[i+1,1] <- paste(format(round(metaresult[1,1], digits = 3), nsmall = 3), 
                                " (", format(round(metaresult[1,2], digits = 3), nsmall = 3), ")", 
                                ifelse(metaresult[1,6] < 0.05, "*", ""), sep = "")
    }
  }
  
  #Fill in the 'TOTAL' column for outcome
  for (i in 1:length(outcomes)){
    subdata <- subset(data, data$outcome == outcomes[i])
    
    #if there are no estimates, then fill the table in as an NA
    if (nrow(subdata) == 0){
      estimates[1,i+1] <- NA
      
    } else if (nrow(subdata) == 1){
      estimates[1, i+1] <- paste(round(subdata$d.directional, 2), " (", round(subdata$sd.d, 2), ")", sep = "")
    } else {
      metaresult <- meta3L(y=subdata$d.directional, v=subdata$var.d, cluster=subdata$article.id)
      metaresult <- summary(metaresult)$coefficients
      
      #store estimate, se, and significance in the estimates table
      estimates[1,i+1] <- paste(format(round(metaresult[1,1], digits = 3), nsmall = 3), 
                                " (", format(round(metaresult[1,2], digits = 3), nsmall = 3), ")", 
                                ifelse(metaresult[1,6] < 0.05, "*", ""), sep = "")
    }
  }
  
  #Fill in the overall total
  metaresult <- meta3L(y=data$d.directional, v=data$var.d, cluster=data$article.id)
  metaresult <- summary(metaresult)$coefficients
  
  #store estimate, se, and significance in the estimates table
  estimates[1,1] <- paste(format(round(metaresult[1,1], digits = 3), nsmall = 3), 
                          " (", format(round(metaresult[1,2], digits = 3), nsmall = 3), ")", 
                          ifelse(metaresult[1,6] < 0.05, "*", ""), sep = "")
  
  return(estimates)
}

# C. Function that extracts relevant estimates and SE's from the results table for coefficient plots
graphable_estimates <- function(estimates, exclude_mechs = F, row = NULL, col = NULL){
  
  ## Logic: the fill_metatable3() function gets us an estimate, SE, and significance decision.
  #         To plot these, we will need to separate things. Let's extract the relevant info:
  
  #Do we want rows or columns?
  if (is.null(row) == T){
    df <- estimates %>% dplyr::select(names(estimates)[col])
  } else if (is.null(col) == T){
    df <- data.frame(t(estimates[row,]))
  }
  names(df)[1] <- "Full"
  
  for(i in 1:nrow(df)){
    myEst <- str_split(df[i, 1], " ")[[1]][1] # get estimate
    mySE <- str_split(df[i, 1], " ")[[1]][2] # get SE
    if(grepl("*", mySE, fixed = T) == T){mySE <- str_sub(mySE, 2, (nchar(mySE)-2))} #remove star
    else {mySE <- str_sub(mySE, 2, (nchar(mySE) - 1))}
    
    df$estimate[i] <- as.numeric(myEst)
    df$se[i] <- as.numeric(mySE)
  }
  
  ## label the estimates that we are getting
  df$var <- row.names(df)
  
  if (is.null(row) == T){
    df$var <- as.factor(df$var)
    df$var <- factor(df$var, levels = c("All victimhood", "General", "Exclusive", "Competitive", "Inclusive"))
    
  } else if (is.null(col) == T){
    df$var <- as.factor(df$var)
    
    if (exclude_mechs == F){
      df$var <- case_match(df$var,
                           "All outcomes" ~ "All\n outcomes", 
                           "Hawkishness" ~ "Hawkishness", 
                           "Reconciliation" ~ "Reconciliation", 
                           "Outgroup Exclusion" ~ "Outgroup \n Exclusion",
                           "Ingroup Attachment" ~ "Ingroup \n Attachment", 
                           "Perspectives" ~ "Perspectives", 
                           "Emotions" ~ "Emotions")
      df$var <- factor(df$var, levels = c("All\n outcomes", "Hawkishness", "Reconciliation", "Outgroup \n Exclusion",
                                          "Ingroup \n Attachment", "Perspectives", "Emotions"))
    } else {
      df$var <- case_match(df$var,
                           "All outcomes" ~ "All\n outcomes", 
                           "Hawkishness" ~ "Hawkishness", 
                           "Reconciliation" ~ "Reconciliation", 
                           "Outgroup Exclusion" ~ "Outgroup \n Exclusion",
                           "Ingroup Attachment" ~ "Ingroup \n Attachment")
      df$var <- factor(df$var, levels = c("All\n outcomes", "Hawkishness", "Reconciliation", "Outgroup \n Exclusion",
                                          "Ingroup \n Attachment"))
    }

  }
  return(df)
}

# D. Function that gives number of estimates and manuscripts in each combination (CVB, outcome) category
get_n <- function(data, exclude_mechs = F, forgraph = F){
  
  ## Logic: we want to know how many manuscripts and how many estimates went into a meta-estimate
  #         Let's extract that information for descriptive statistics and for in-plot text
  
  if (exclude_mechs == F){
    #1. get number of estimates
    ntab <- data.frame(unclass(table(data$treatment, data$outcome)))
    ntab$iv.total <- rowSums(ntab)
    ntab <- ntab[,c(3, 6, 5, 4, 1, 2, 7)]
    ntab <- rbind(ntab, dv.total = colSums(ntab))

    #2. get number of manuscripts
    unique_est <- unique(data[,c('article.id', 'treatment', 'outcome')])
    mtab <- data.frame(unclass(table(unique_est$treatment, unique_est$outcome)))
    mtab <- rbind(mtab, dv.total = table(unique(data[,c('article.id', 'outcome')])$outcome))
    mtab$iv.total <- c(table(unique(data[,c('article.id', 'treatment')])$treatment), NA)
    mtab <- mtab[,c(3, 6, 5, 4, 1, 2, 7)]
    mtab[5,7] <- length(unique(unique_est$article.id))
    
    #3. create combined table
    combotab <- ntab
    for (i in 1:5){for (j in 1:7){
      combotab[i,j] <- paste(ntab[i,j], " (", mtab[i,j], ")", sep = "")
    } }
  } else{
    
    #1. get number of estimates
    ntab <- data.frame(unclass(table(data$treatment, data$outcome)))
    ntab$iv.total <- rowSums(ntab)
    ntab <- ntab[,c(1, 4, 3, 2, 5)]
    ntab <- rbind(ntab, dv.total = colSums(ntab))
    
    #2. get number of manuscripts
    unique_est <- unique(data[,c('article.id', 'treatment', 'outcome')])
    mtab <- data.frame(unclass(table(unique_est$treatment, unique_est$outcome)))
    mtab <- rbind(mtab, dv.total = table(unique(data[,c('article.id', 'outcome')])$outcome))
    mtab$iv.total <- c(table(unique(data[,c('article.id', 'treatment')])$treatment), NA)
    mtab <- mtab[,c(1, 4, 3, 2, 5)]
    mtab[5,5] <- length(unique(unique_est$article.id))
    
    #3. create combined table
    combotab <- ntab
    for (i in 1:5){for (j in 1:5){
      combotab[i,j] <- paste(ntab[i,j], " (", mtab[i,j], ")", sep = "")
    } }
    
  }
  
  row.names(ntab) <- c("Generic", "Competitive", "Exclusive", "Inclusive", "Total")
  row.names(mtab) <- c("Generic", "Competitive", "Exclusive", "Inclusive", "Total")
  row.names(combotab) <- c("Generic", "Competitive", "Exclusive", "Inclusive", "Total")
  
  if (exclude_mechs == T){
    if (forgraph == T){ #Special formatting for full_plots graph
      combotab <- combotab[c(5,1:4),]
      combotab <- combotab[, c(5,1:4)]
      }
  } else {
    if (forgraph == T){ 
      combotab <- combotab[c(5, 1:4),]
      combotab <- combotab[, c(7, 1:6)]
      }
  }
  
  return(list(ntab = ntab, mtab = mtab, combotab = combotab))
}

# E. Function that makes the coefficient plot for each individual estimate 
full_plots <- function(metatable, df = data, exclude_mechs = F, range = 1.7, jump = 0.5, Bonf = TRUE, only_mech = F ){
  
  #get meta estimates
  main <- metatable
  if (exclude_mechs == F) {ncolumn = 7
  } else {ncolumn = 5}
  
  
  #make graphable table from these estimates
  ind_coefs <- rbind(
    cbind(graphable_estimates(main, exclude_mechs = exclude_mechs, row = 1), "Victimhood Belief" = rep("All Types", ncolumn)),
    cbind(graphable_estimates(main, exclude_mechs = exclude_mechs, row = 2), "Victimhood Belief" = rep("Generic", ncolumn)),
    cbind(graphable_estimates(main, exclude_mechs = exclude_mechs, row = 3), "Victimhood Belief" = rep("Competitive", ncolumn)),
    cbind(graphable_estimates(main, exclude_mechs = exclude_mechs, row = 4), "Victimhood Belief" = rep("Exclusive", ncolumn)),
    cbind(graphable_estimates(main, exclude_mechs = exclude_mechs, row = 5), "Victimhood Belief" = rep("Inclusive", ncolumn))
  )
  
  
  ### What kind of multiple-hypothesis correction are we doing (can choose Bonferroni (Bonf == T) or Benjamini-Hochberg)?
  
  if(Bonf == TRUE){
          ### Bonferroni Correction
    
    #make adjustable SE column (will change for where we only make 1, 4, or 6 comparisons)
    ind_coefs$special_se <- ind_coefs$se
    
    if (exclude_mechs == F){ #If we are including cognitive frames and emotions
      ind_coefs[1,6] <- ind_coefs[1,6]*(1.645/1.645)  #Full estimate [1 correction]
      ind_coefs[2:7,6] <- ind_coefs[2:7,6]*(2.394/1.645) #All victimhood [6 corrections]
      ind_coefs[c(8, 15, 22, 29),6] <- ind_coefs[c(8, 15, 22, 29),6]*(2.2414/1.645) #All Outcomes [4 corrections]
      ind_coefs[-c(1:8, 15, 22, 29), 6] <- ind_coefs[-c(1:8, 15, 22, 29), 6]*(2.865/1.645) #Individual Estimates [24 corrections]
    } else {
      ind_coefs[1,6] <- ind_coefs[1,6]*(1.645/1.645)  #Full estimate [1 correction]
      ind_coefs[2:5,6] <- ind_coefs[2:5,6]*(2.2414/1.645) #All victimhood [4 correction]
      ind_coefs[c(6, 11, 16, 21),6] <- ind_coefs[c(6, 11, 16, 21),6]*(2.2414/1.645) #All Outcomes [4 correction]
      ind_coefs[-c(1:6, 11, 16, 21), 6] <- ind_coefs[-c(1:6, 11, 16, 21), 6]*(2.7344/1.645) #Individual Estimates [16 corrections]
    }
    
    if (only_mech == T){
      ind_coefs$special_se <- ind_coefs$se
      ind_coefs[6:7, 6] <- ind_coefs[6:7, 6]*(1.996/1.645) #All victimhood [2 corrections]
      ind_coefs[c(13, 14, 20, 21, 27, 28, 34, 35), 6] <- ind_coefs[c(13, 14, 20, 21, 27, 28, 34, 35), 6]*(2.498/1.645) #ind est [8 corr]
    }

    
  } else{
        ### Benjamini-Hochberg Correction
    
    #make adjustable SE column (will change for where we only make 1, 4, or 6 comparisons)
    ind_coefs$special_se <- ind_coefs$se
    
    if (exclude_mechs == F){
      #Full estimate
      ind_coefs[1,6] <- ind_coefs[1,6]*(1.645/1.645) 
      
      #All victimhood
      ranking = rank(pnorm(-abs(ind_coefs[2:7, 2]/ind_coefs[2:7, 6]), 0, 1), ties.method = 'first')
      ind_coefs[2:7,6] <- ind_coefs[2:7,6]*(abs(qnorm(ranking*0.05/length(ranking)))/1.645) 
      
      #All Outcomes
      ranking = rank(pnorm(-abs(ind_coefs[c(8, 15, 22, 29), 2]/ind_coefs[c(8, 15, 22, 29), 6]), 0, 1), ties.method = 'first')
      ind_coefs[c(8, 15, 22, 29), 6] <- ind_coefs[c(8, 15, 22, 29),6]*(abs(qnorm(ranking*0.05/length(ranking)))/1.645)
      
      
      #Individual estimates
      ranking = rank(pnorm(-abs(ind_coefs[-c(1:8, 15, 22, 29), 2]/ind_coefs[-c(1:8, 15, 22, 29), 6]), 0, 1), ties.method = 'first')
      ind_coefs[-c(1:8, 15, 22, 29), 6] <- ind_coefs[-c(1:8, 15, 22, 29), 6]*(abs(qnorm(ranking*0.05/length(ranking)))/1.645) 
    } else {
      #Full estimate
      ind_coefs[1,6] <- ind_coefs[1,6]*(1.645/1.645) 
      
      #All victimhood
      ranking = rank(pnorm(-abs(ind_coefs[2:5, 2]/ind_coefs[2:5, 6]), 0, 1), ties.method = 'first')
      ind_coefs[2:5,6] <- ind_coefs[2:5,6]*(abs(qnorm(ranking*0.05/length(ranking)))/1.645) 
      
      #All Outcomes
      ranking = rank(pnorm(-abs(ind_coefs[c(6, 11, 16, 21), 2]/ind_coefs[c(6, 11, 16, 21), 6]), 0, 1), ties.method = 'first')
      ind_coefs[c(6, 11, 16, 21), 6] <- ind_coefs[c(6, 11, 16, 21),6]*(abs(qnorm(ranking*0.05/length(ranking)))/1.645)
      
      
      #Individual estimates
      ranking = rank(pnorm(-abs(ind_coefs[-c(1:6, 11, 16, 21), 2]/ind_coefs[-c(1:6, 11, 16, 21), 6]), 0, 1), ties.method = 'first')
      ind_coefs[-c(1:6, 11, 16, 21), 6] <- ind_coefs[-c(1:6, 11, 16, 21), 6]*(abs(qnorm(ranking*0.05/length(ranking)))/1.645) 
    }
    

  }

  
  
  #Change variable name on Reconciliation to 'Opposition to \n Reconciliation
  ind_coefs$var <- case_match(ind_coefs$var, "Reconciliation" ~ "Opposition to \n Reconciliation", .default = ind_coefs$var)
  ind_coefs$var <- case_match(ind_coefs$var, "All\n outcomes" ~ "All \n Outcomes", .default = ind_coefs$var)
  
  if (exclude_mechs == F){
    ind_coefs$var <- factor(ind_coefs$var, levels = c("All \n Outcomes", "Hawkishness", "Opposition to \n Reconciliation", 
                                                      "Outgroup \n Exclusion","Ingroup \n Attachment", "Perspectives", "Emotions"))
  } else {
    ind_coefs$var <- factor(ind_coefs$var, levels = c("All \n Outcomes", "Hawkishness", "Opposition to \n Reconciliation", 
                                                      "Outgroup \n Exclusion","Ingroup \n Attachment"))
  }

  
  #Add number of estimates and manuscripts
  ind_coefs$n <- rep(NA, nrow(ind_coefs))
  for (i in 1:5){
    for (j in 1:ncolumn){
      ind_coefs$n[((i-1)*ncolumn)+j] <- paste("[", get_n(df, exclude_mechs = exclude_mechs, forgraph = T)$combotab[i,j],"]", sep ="")
    }
  }
  ind_coefs$char <- nchar(ind_coefs$n) - 2 #characters are for sizing/placing numbers
  
  #make the victimhood type factor the right order
  ind_coefs$`Victimhood Belief` <- factor(ind_coefs$`Victimhood Belief`,
                                        levels = c("All Types", "Generic", "Competitive", "Exclusive", "Inclusive"))
  
  if (only_mech == T){
    ind_coefs <- ind_coefs[c(6, 7, 13, 14, 20, 21, 27, 28, 34, 35),] #subset to just emotions and cog. frames
  }
  
  #make graph
  final_plot <- ggplot(data = ind_coefs, aes(label = n, shape = `Victimhood Belief`, color = `Victimhood Belief`, y = var))+
    geom_pointrange(aes(x = estimate, xmax = estimate + 1.645*special_se, xmin = estimate - 1.645*special_se),linewidth = 0.7,
                    position = position_dodge(width = -1/2)) +
    geom_pointrange(aes(x = estimate, xmax = estimate + 1.645*se, xmin = estimate - 1.645*se), size = 0.9, linewidth = 1.2,
                    position = position_dodge(width = -1/2)) +
    
    scale_shape_manual(values = c(15, 16, 17, 4, 3)) +
    scale_color_manual(values = c("black", "red", "blue", "darkgreen", "darkorange")) +
    
    geom_text(aes(x = ifelse(`Victimhood Belief` == "Inclusive", 
                             estimate - 3*se - 0.07*((char-1)/3), estimate + 3*se + 0.07*((char-1)/3))),
              fontface = "bold", position = position_dodge(width = -1/2), size = 4, show_guide = F) +
    
    geom_vline(xintercept = 0, linetype = "dashed")+
    scale_y_discrete(limits = rev) +
    coord_cartesian(xlim = c(-range, range), clip = 'off') +
    scale_x_continuous(breaks = seq(round(-range), round(range), jump)) +
    labs(x = "Effect size meta-estimate (90% CI)") +
    theme_bw() +
    theme(axis.title.y = element_blank(),
          axis.text = element_text(size = 14),
          legend.text = element_text(size = 14),
          axis.title.x = element_text(size = 14),
          legend.title = element_text(size = 14))
  return(final_plot)
}

# F. Function that makes coefficient plot to compare across a binary split (for example: experimental vs. observational data)
comparison_plot <- function(data1, data2, df1, df2, name1 = "Experimental", name2 = "Observational",
                         exclude_mechs = T, range = 1, jump = 0.5){
  
  #create datatable to hold information
  tab <- rbind(
    cbind(graphable_estimates(data1, col = 1), type = rep(name1, 5), 
          n = get_n(df1, exclude_mechs = exclude_mechs, forgraph = T)$combotab[, 1]),
    
    cbind(graphable_estimates(data2, col = 1), type = rep(name2, 5), 
          n = get_n(df2, exclude_mechs = exclude_mechs, forgraph = T)$combotab[,1]))
  
  tab$n <- paste("[", tab$n, "]", sep = "")
  tab$char <- nchar(tab$n) - 2 #characters are for sizing/placing numbers
  
  tab$var <- case_match(tab$var, "General" ~ "Generic", "All victimhood" ~ "All victimhood\nbeliefs",  .default = tab$var)
  tab$var <- factor(tab$var, levels = c("All victimhood\nbeliefs", "Generic", "Competitive", "Exclusive", "Inclusive"))
  
  mylabs <- c( ' Number of estimates, manuscripts\n in brackets [n (m)]' )
  
  #plot results
  myplot <- ggplot(data = tab, aes(y = var, shape = type, color = type, label = n))+
    
    geom_pointrange(aes(x = estimate, 
                        xmax = estimate + 1.645*se, xmin = estimate - 1.645*se, 
                        y = var),lwd = 1.2, size = 0.9,
                    position = position_dodge(width = -1/2)) +
    geom_text(aes(x = ifelse(var == "Inclusive", 
                             estimate - 2.05*se - 0.05*((char-1)/3), estimate + 2.05*se + 0.05*((char-1)/3))),
              fontface = "bold", position = position_dodge(width = -1/2), size = 4, show_guide = F) +
    
    scale_shape_manual(values = c(15, 17)) +
    scale_color_manual(values = c("darkblue", "darkred")) +
    geom_vline(xintercept = 0, linetype = "dashed")+
  
    scale_y_discrete(limits = rev) +
    scale_x_continuous(limits = c(-range, range), breaks = seq(round(-range, 1), round(range, 1), jump)) +
    labs(x = "Effect size meta-estimate (90% CI)") +
    theme_bw() +
    theme(axis.title.y = element_blank(),
          axis.title.x = element_text(size = 14),
          axis.text = element_text(size = 14),
          legend.text = element_text(size = 14),
          legend.title = element_blank())
  return(myplot)
}

# G. Write a function that calculates coefficients when examining sub_samples
subset_estimates <- function(df, variable){
  
  data <- subset(df, is.na(df[,variable]) == F)
  row.names(data) <- seq(1, nrow(data), 1) #re-index
  
  #1. Non-Inclusive Victimhood
  data1 <- subset(data, data$treatment != "Inclusive victimhood")
  data1$relevant_var1 <- data1[,variable] #save relevant variable
  results1 <- data.frame(setting = sort(unique(data1$relevant_var1)),
                         estimate = rep(NA, length(unique(data1$relevant_var1))),
                         se = rep(NA, length(unique(data1$relevant_var1))),
                         n = paste("[", table(data1$relevant_var1), " (", 
                                   table(unique(data.frame(v1 = data1$relevant_var1, v2 = data1$article.id))$v1),
                                   ")]", sep = ""),
                         type = rep("Non-Inclusive victimhood", length(unique(data1$relevant_var1))))
  
  #1.1 Add estimates and Standard Errors
  for (i in 1:nrow(results1)){
    df <- create_weights(subset(data1, data1$relevant_var1 == results1$setting[i]))
    
    if(nrow(df) == 1){ #if only one row, cannot run meta3L. Just take estimate and sd
      results1[i, 2:3] <- c(df$d.directional, df$sd.d)
    } else{
      metaresult <- meta3L(y = df$d.directional, v = df$var.d, cluster = df$article.id)
      metaresult <- summary(metaresult)$coefficients
      results1[i, 2:3] <- c(metaresult[1,1], metaresult[1,2])}
  }
  
  
  #2. Inclusive Victimhood
  data2 <- subset(data, data$treatment == "Inclusive victimhood")
  data2$relevant_var2 <- data2[,variable] #save setting variable
  
  results2 <- data.frame(setting = sort(unique(data2$relevant_var2)),
                         estimate = rep(NA, length(unique(data2$relevant_var2))),
                         se = rep(NA, length(unique(data2$relevant_var2))),
                         n = paste("[", table(data2$relevant_var2), " (", 
                                   table(unique(data.frame(v1 = data2$relevant_var2, v2 = data2$article.id))$v1),
                                   ")]", sep = ""),
                         type = rep("Inclusive victimhood", length(unique(data2$relevant_var2))))
  
  #2.1 Add estimates and Standard Errors
  for (i in 1:nrow(results2)){
    df <- create_weights(subset(data2, data2$relevant_var2 == results2$setting[i]))
    
    if(nrow(df) == 1){ #if only one row, cannot run meta3L. Just take estimate and sd
      results2[i, 2:3] <- c(df$d.directional, df$sd.d)
    } else{
      metaresult <- meta3L(y = df$d.directional, v = df$var.d, cluster = df$article.id)
      metaresult <- summary(metaresult)$coefficients
      results2[i, 2:3] <- c(metaresult[1,1], metaresult[1,2])}
  }

  #3. Bind into a single dataframe and return
  results <- rbind(results1, results2)
  results$char <- nchar(results$n) - 2 #characters are for sizing/placing numbers
  
  return(results)
}

# H. Write a function that plots subset shift estimates
subset_plots <- function(subset_data, range = 2.2, jump = 1, separate = F, sample_size = F, ind_plots = F){
  
  #Get actual meta-estimates for plot subtitle
  incdata <- subset(main, main$treatment == "Inclusive victimhood")
  INCmeta <- meta3L(y = incdata$d.directional, v = incdata$var.d, cluster = incdata$article.id)
  inc.est <- round(summary(INCmeta)$coefficients[1,1], 2)
  
  exdata <- subset(main, main$treatment != "Inclusive victimhood")
  EXmeta <- meta3L(y = exdata$d.directional, v = exdata$var.d, cluster = exdata$article.id)
  ex.est <- round(summary(EXmeta)$coefficients[1,1], 2)
  
  #If we want separate graphs...
  if (separate == T){
    non_uni <- ggplot(data = subset(subset_data, subset_data$type != "Inclusive victimhood")) +
      geom_pointrange(aes(x = estimate, y = reorder(setting, estimate), 
                          xmin = estimate - 1.645*se,
                          xmax = estimate + 1.645*se), alpha = 1, color = "darkred", shape = 17, lwd = 1.2) +
      geom_vline(xintercept = 0, linetype = "dashed")+
      geom_vline(xintercept = ex.est, color = "darkred", alpha = 0.5, lwd = 1) +
      geom_text(aes(x = estimate + 2.2*se + 0.1*((char-1)/3), y = setting, label = n),
                color = "darkred", fontface = 'bold', size = 4, show_guide = F)+
      
      labs(x = "Effect size meta-estimate (90% CI)",
           subtitle = paste("Meta estimate, d = ", ex.est, sep = "")) +
      scale_x_continuous(limits = c(-range, range), breaks = seq(-round(range), round(range), jump)) +
      theme_bw() +
      theme(axis.title.y = element_blank(),
            axis.ticks.y = element_blank(),
            panel.grid.major.y = element_blank(),
            axis.text.y = element_text(size = 12),
            plot.subtitle = element_text(size = 14),
            axis.title.x = element_text(size = 14))
    
    
    uni <- ggplot(data = subset(subset_data, subset_data$type == "Inclusive victimhood")) +
      geom_pointrange(aes(x = estimate, y = reorder(setting, estimate), 
                          xmin = estimate - 1.645*se,
                          xmax = estimate + 1.645*se), alpha = 1, shape = 15, color = "darkblue", lwd = 1.2) +
      geom_vline(xintercept = 0, linetype = "dashed")+
      geom_vline(xintercept = inc.est, color = "darkblue", alpha = 0.5, lwd = 1) +
      geom_text(aes(x = estimate - 2.2*se - 0.1*((char-1)/3), y = setting, label = n),
                color = "darkblue", fontface = 'bold', size = 4, show_guide = F)+
      
      labs(x = "Effect size meta-estimate (90% CI)",
           subtitle = paste("Meta estimate, d = ", inc.est, sep = "")) +
      scale_x_continuous(limits = c(-range, range), breaks = seq(-round(range), round(range), jump)) +
      theme_bw() +
      theme(axis.title.y = element_blank(),
            axis.ticks.y = element_blank(),
            panel.grid.major.y = element_blank(),
            axis.text.y = element_text(size = 12),
            plot.subtitle = element_text(size = 14),
            axis.title.x = element_text(size = 14))
    
    if(ind_plots == T){
      myplot <- list(uni = uni, non_uni = non_uni)
    }else{
      myplot <- uni + non_uni
    }
    
  } else {
    
    if (sample_size == T){
      subset_data$setting <- factor(subset_data$setting, 
                                    levels = c("[25,110]", "(110,225]", "(225,461]", "(461, 29169]"))}
    
    subset_data$type <- ifelse(subset_data$type == "Inclusive victimhood", "\nInclusive\nvictimhood\nbeliefs\n", 
                               "Non-Inclusive\nvictimhood\nbeliefs")
    
    #(If we want one graph with offset estimates)
    myplot <- ggplot(data = subset_data, aes(x = estimate, y = setting, label = n, shape = type, color = type)) +
      geom_pointrange(aes(
        xmin = estimate - 1.645*se,
        xmax = estimate + 1.645*se), alpha = 1, size = 0.9, lwd = 1.2,
        position = position_dodge(width = -1/2)) +
      geom_vline(xintercept = 0, linetype = "dashed")+
      geom_vline(xintercept = ex.est, color = "darkred", alpha = 0.5, lwd = 1) +
      geom_vline(xintercept = inc.est,  color = "darkblue", alpha = 0.5, lwd = 1) +
      
      labs(x = "Effect size meta-estimate (90% CI)",
           subtitle = paste("Meta estimates as solid lines: d = ", inc.est, " (inclusive), d = ", ex.est, " (non-inclusive)", sep = "")) +
      scale_shape_manual(values = c(15, 17)) + 
      scale_color_manual(values = c("darkblue", "darkred"))+
      
      scale_y_discrete(limits = rev)+
      scale_x_continuous(limits = c(-range, range), breaks = seq(-round(range), round(range), jump)) +
      theme_bw() +
      theme(axis.title.y = element_blank(),
            axis.ticks.y = element_blank(),
            panel.grid.major.y = element_blank(),
            axis.text.y = element_text(size = 12),
            plot.subtitle = element_text(size = 14),
            axis.title.x = element_text(size = 14),
            legend.text = element_text(size = 14),
            legend.title = element_blank()) #,
            #legend.position = "bottom")
    
    if (sample_size ==F){
      myplot <- myplot + 
        geom_text(aes(x = ifelse(type == "Inclusive victimhood", 
                                 estimate - 2.2*se - 0.1*((char-1)/3), estimate + 2.2*se + 0.1*((char-1)/3))),
                  fontface = "bold", position = position_dodge(width = -1/2), size = 4, show_guide = F)
    }
  }
  
  return(myplot)
  
}





############## APPENDIX ONLY ####
# A1. Create function that adds weighs to data  (can be by estimates/manuscript or by inverse variance weighting)
create_weights <- function(data, ivw = F, cap = 10000, stage2ivw = F){
  
  #remove existing weights column if there
  if ("weights" %in% colnames(data)){ data <- data[, ! colnames(data) %in% c("weights")]}
  
  #Weights by number of estimates/manuscript
  if(ivw==F){
    weight_data <- data %>%
      group_by(Title) %>%
      summarise(weights = 1/n()) %>%
      as.data.frame()
    
    data <- merge(data, weight_data, by = "Title")
    
  } else{ #Weights from inverse variance
    if (stage2ivw == F){data$weights <- 1/data$var.d #take regular variance
    } else { data$weights <- 1/data$v_ind} #take variance from two stage calculation
    
    for(i in 1:nrow(data)){
      data$weights[i] <- ifelse(data$weights[i]>cap, cap, data$weights[i])
    }
  }
  return(data)
}

# A2. Write function that fills table for meta-estimates (d and se) -- 2-stage simple REML
fillMeta_simple <- function(data, exclude_mechs = F, weights = T, ivw = F, cap = 10000, model = 'REML'){
  
  #### A. create dataframe to store results
  
  if (exclude_mechs == F){ # if we are including emotions and cognitive frames...
    estimates <- as.data.frame(matrix(nrow = 5, ncol = 7))
    row.names(estimates) <- c("All victimhood", "General", "Competitive", "Exclusive", "Inclusive")
    names(estimates) <- c("All outcomes","Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment", "Perspectives", "Emotions")
    
    #specify the victimhood types and outcomes
    victimhood <- c("Collective victimhood", "Competitive victimhood", "Exclusive victimhood", "Inclusive victimhood")
    outcomes <- c("Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment", "Cognitive Frame", "Emotions")
  } else {
    estimates <- as.data.frame(matrix(nrow = 5, ncol = 5))
    row.names(estimates) <- c("All victimhood", "General", "Competitive", "Exclusive", "Inclusive")
    names(estimates) <- c("All outcomes","Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment")
    
    #specify the victimhood types and outcomes
    victimhood <- c("Collective victimhood", "Competitive victimhood", "Exclusive victimhood", "Inclusive victimhood")
    outcomes <- c("Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment")
  }
  
  
  #### B. Make table 4 x 6 (5x7 with totals) table [or 4 x 4 (5x5 with totals)]
  for (i in 1:length(victimhood)){
    for (j in 1:length(outcomes)){
      #subset by victimhood category
      subdata <- subset(data, data$treatment == victimhood[i])
      
      #subset by outcome category
      subdata <- subset(subdata, subdata$outcome == outcomes[j])
      
      #if there are no estimates, then fill the table in as an NA
      if (nrow(subdata) == 0){
        estimates[i+1,j+1] <- NA
      } else {
        #create weights for this particular combination
        subdata <- create_weights(subdata, ivw = ivw, cap = cap)
        
        #run meta-analysis
        if (weights == T){
          metaresult <- rma(yi = subdata$d.directional, vi = subdata$var.d,
                            weights = subdata$weights, method = model)
        } else{
          metaresult <- rma(yi = subdata$d.directional, vi = subdata$var.d, method = model)
        }
        
        #store estimate, se, and significance in the estimates table
        estimates[i+1,j+1] <- paste(format(round(metaresult$b[,1], digits = 3), nsmall = 3), 
                                    " (", format(round(metaresult$se, digits = 3), nsmall = 3), ")", 
                                    ifelse(metaresult$pval < 0.05, "*", ""), sep = "") 
      }
    }
  }
  
  
  #Fill in the 'TOTAL' column for victimood cases
  for (i in 1:length(victimhood)){
    subdata <- subset(data, data$treatment == victimhood[i])
    
    #if there are no estimates, then fill the table in as an NA
    if (nrow(subdata) == 0){
      estimates[i+1,1] <- NA
      
    } else {
      subdata <- create_weights(subdata, ivw = ivw, cap = cap)
      
      if (weights == T){
        metaresult <- rma(yi = subdata$d.directional, vi = subdata$var.d,
                          weights = subdata$weights, method = model)
      } else{
        metaresult <- rma(yi = subdata$d.directional, vi = subdata$var.d, method = model)
      }
      
      estimates[i+1,1] <- paste(format(round(metaresult$b[,1], digits = 3), nsmall = 3), 
                                " (", format(round(metaresult$se, digits = 3), nsmall = 3), ")", 
                                ifelse(metaresult$pval < 0.05, "*", ""), sep = "") 
    }
  }
  
  #Fill in the 'TOTAL' column for victimood cases
  for (i in 1:length(outcomes)){
    subdata <- subset(data, data$outcome == outcomes[i])
    
    #if there are no estimates, then fill the table in as an NA
    if (nrow(subdata) == 0){
      estimates[1,i+1] <- NA
      
    } else {
      subdata <- create_weights(subdata, ivw = ivw, cap = cap)
      
      if (weights == T){
        metaresult <- rma(yi = subdata$d.directional, vi = subdata$var.d,
                          weights = subdata$weights, method = model)
      } else{
        metaresult <- rma(yi = subdata$d.directional, vi = subdata$var.d, method = model)
      }
      
      estimates[1,i+1] <- paste(format(round(metaresult$b[,1], digits = 3), nsmall = 3), 
                                " (", format(round(metaresult$se, digits = 3), nsmall = 3), ")", 
                                ifelse(metaresult$pval < 0.05, "*", ""), sep = "")
    }
  }
  
  #Fill in the overall total
  data <- create_weights(data, ivw = ivw, cap = cap)
  if (weights == T){
    metaresult <- rma(yi = data$d.directional, vi = data$var.d,
                      weights = data$weights, method = model)
  } else{
    metaresult <- rma(yi = data$d.directional, vi = data$var.d, method = model)
  }
  
  estimates[1,1] <- paste(format(round(metaresult$b[,1], digits = 3), nsmall = 3), 
                          " (", format(round(metaresult$se, digits = 3), nsmall = 3), ")", 
                          ifelse(metaresult$pval < 0.05, "*", ""), sep = "")
  
  return(estimates)
}

# J. Write function that fills in the 6x7 table using a 2-stage meta analysis process
fillMeta_twostage <- function(data, exclude_mechs = F, iv_weights = F, cap = 10000, model = 'REML'){
  
  #### FIRST STEP: Get d and variance estimates for each UNIQUE article (first stage)
  stage2 <- cbind(data, matrix(data = NA, nrow = nrow(data), ncol = 8))
  names(stage2)[(ncol(data) + 1):(ncol(data) + 8)] <- c("d_all", "v_all", "d_IV", "v_IV", "d_DV", "v_DV", "d_ind", "v_ind")
  
  
  #for the article.id grouping (all)
  for (i in 1:length(unique(stage2$article.id))){
    df <- subset(stage2, stage2$article.id == unique(stage2$article.id)[i])
    
    metaresult <- rma(yi = df$d.directional, vi = df$var.d, method = "REML")
    stage2[c(row.names(df)), c('d_all')] <- as.numeric(metaresult$b)
    stage2[c(row.names(df)), c('v_all')] <- (as.numeric(metaresult$se))^2
    
    # For article ID grouping + Victimhood Category
    for (j in 1:length(unique(df$treatment))){
      df_IV <- subset(df, df$treatment == unique(df$treatment)[j])
      metaresult <- rma(yi = df_IV$d.directional, vi = df_IV$var.d, method = "REML")
      stage2[c(row.names(df_IV)), c('d_IV')] <- as.numeric(metaresult$b)
      stage2[c(row.names(df_IV)), c('v_IV')] <- (as.numeric(metaresult$se))^2
    }
    
    # For article ID grouping + Outcome Category
    for (k in 1:length(unique(df$outcome))){
      df_DV <- subset(df, df$outcome == unique(df$outcome)[k])
      metaresult <- rma(yi = df_DV$d.directional, vi = df_DV$var.d, method = "REML")
      stage2[c(row.names(df_DV)), c('d_DV')] <- as.numeric(metaresult$b)
      stage2[c(row.names(df_DV)), c('v_DV')] <- (as.numeric(metaresult$se))^2
      
      #For each individual article + Victimhood type + Outcome Category
      for (x in 1: length(unique(df_DV$treatment))){
        df_ind <- subset(df_DV, df_DV$treatment == unique(df_DV$treatment)[x])
        metaresult <- rma(yi = df_ind$d.directional, vi = df_ind$var.d, method = "REML")
        stage2[c(row.names(df_ind)), c('d_ind')] <- as.numeric(metaresult$b)
        stage2[c(row.names(df_ind)), c('v_ind')] <- (as.numeric(metaresult$se))^2
      }
    }
    
  }
  
  
  #### A. create dataframe to store results
  
  if (exclude_mechs == F){ # if we are including emotions and cognitive frames...
    estimates <- as.data.frame(matrix(nrow = 5, ncol = 7))
    row.names(estimates) <- c("All victimhood", "General", "Competitive", "Exclusive", "Inclusive")
    names(estimates) <- c("All outcomes","Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment", "Perspectives", "Emotions")
    
    #specify the victimhood types and outcomes
    victimhood <- c("Collective victimhood", "Competitive victimhood", "Exclusive victimhood", "Inclusive victimhood")
    outcomes <- c("Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment", "Cognitive Frame", "Emotions")
  } else {
    estimates <- as.data.frame(matrix(nrow = 5, ncol = 5))
    row.names(estimates) <- c("All victimhood", "General", "Competitive", "Exclusive", "Inclusive")
    names(estimates) <- c("All outcomes","Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment")
    
    #specify the victimhood types and outcomes
    victimhood <- c("Collective victimhood", "Competitive victimhood", "Exclusive victimhood", "Inclusive victimhood")
    outcomes <- c("Hawkishness", "Reconciliation", "Outgroup Exclusion", "Ingroup Attachment")
  }
  
  
  #### B. Make table 4 x 6 (5x7 with totals) table [or 4 x 4 (5x5 with totals)]
  for (i in 1:length(victimhood)){
    for (j in 1:length(outcomes)){
      #subset by victimhood category
      subdata <- subset(stage2, stage2$treatment == victimhood[i])
      
      #subset by outcome category
      subdata <- subset(subdata, subdata$outcome == outcomes[j])
      
      subdata <- subdata[match(unique(subdata$article.id), subdata$article.id),]
      
      
      #if there are no estimates, then fill the table in as an NA
      if (nrow(subdata) == 0){
        estimates[i+1,j+1] <- NA
        
        
      } else {
        #create weights for this particular combination
        subdata <- create_weights(subdata, ivw = T, cap = cap, stage2ivw = T)
        
        #run meta-analysis
        if (iv_weights == T){
          metaresult <- rma(yi = subdata$d_ind, vi = subdata$v_ind,
                            weights = subdata$weights, method = model)
        } else{
          metaresult <- rma(yi = subdata$d_ind, vi = subdata$v_ind, method = model)
        }
        
        #store estimate, se, and significance in the estimates table
        estimates[i+1,j+1] <- paste(format(round(metaresult$b[,1], digits = 3), nsmall = 3), 
                                    " (", format(round(metaresult$se, digits = 3), nsmall = 3), ")", 
                                    ifelse(metaresult$pval < 0.05, "*", ""), sep = "") 
      }
    }
  }
  
  
  #Fill in the 'TOTAL' column for victimood cases
  for (i in 1:length(victimhood)){
    subdata <- subset(stage2, stage2$treatment == victimhood[i])
    subdata <- subdata[match(unique(subdata$article.id), subdata$article.id),]
    
    #if there are no estimates, then fill the table in as an NA
    if (nrow(subdata) == 0){
      estimates[i+1,1] <- NA
      
    } else {
      subdata <- create_weights(subdata, ivw = T, cap = cap, stage2ivw = T)
      
      if (iv_weights == T){
        metaresult <- rma(yi = subdata$d_IV, vi = subdata$v_IV,
                          weights = subdata$weights, method = model)
      } else{
        metaresult <- rma(yi = subdata$d_IV, vi = subdata$v_IV, method = model)
      }
      
      
      estimates[i+1,1] <- paste(format(round(metaresult$b[,1], digits = 3), nsmall = 3), 
                                " (", format(round(metaresult$se, digits = 3), nsmall = 3), ")", 
                                ifelse(metaresult$pval < 0.05, "*", ""), sep = "") 
    }
  }
  
  #Fill in the 'TOTAL' column for victimood cases
  for (i in 1:length(outcomes)){
    subdata <- subset(stage2, stage2$outcome == outcomes[i])
    subdata <- subdata[match(unique(subdata$article.id), subdata$article.id),]
    
    #if there are no estimates, then fill the table in as an NA
    if (nrow(subdata) == 0){
      estimates[1,i+1] <- NA
      
    } else {
      subdata <- create_weights(subdata, ivw = T, cap = cap, stage2ivw = T)
      
      if (iv_weights == T){
        metaresult <- rma(yi = subdata$d_DV, vi = subdata$v_DV,
                          weights = subdata$weights, method = model)
      } else{
        metaresult <- rma(yi = subdata$d_DV, vi = subdata$v_DV, method = model)
      }
      
      estimates[1,i+1] <- paste(format(round(metaresult$b[,1], digits = 3), nsmall = 3), 
                                " (", format(round(metaresult$se, digits = 3), nsmall = 3), ")", 
                                ifelse(metaresult$pval < 0.05, "*", ""), sep = "")
    }
  }
  
  #Fill in the overall total
  subdata <- stage2[match(unique(stage2$article.id), stage2$article.id),]
  
  subdata <- create_weights(subdata, ivw = T, cap = cap, stage2ivw = T)
  if (iv_weights == T){
    metaresult <- rma(yi = subdata$d_all, vi = subdata$v_all,
                      weights = subdata$weights, method = model)
  } else{
    metaresult <- rma(yi = subdata$d_all, vi = subdata$v_all, method = model)
  }
  
  estimates[1,1] <- paste(format(round(metaresult$b[,1], digits = 3), nsmall = 3), 
                          " (", format(round(metaresult$se, digits = 3), nsmall = 3), ")", 
                          ifelse(metaresult$pval < 0.05, "*", ""), sep = "")
  
  return(estimates)
}

